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PARAMETER CHOICE STRATEGIES FOR LEAST-SQUARES 
APPROXIMATION OF NOISY SMOOTH FUNCTIONS ON THE 

SPHERE 

S. V. PEREVERZYEV*, I. H. SLOANt, AND P. TKACHENKO* 

Abstract. We consider a polynomial reconstruction of smooth functions from their noisy values 
at discrete nodes on the unit sphere by a variant of the regularized least-squares method of An et 
al., SIAM J. Numer. Anal. 50 (2012), 1513-1534. As nodes we use the points of a positive-weight 
cubature formula that is exact for all spherical polynomials of degree up to 2 M, where M is the 
degree of the reconstructing polynomial. We first obtain a reconstruction error bound in terms of 
the regularization parameter and the penalization parameters in the regularization operator. Then we 
discuss a priori and a posteriori strategies for choosing these parameters. Finally, we give numerical 
examples illustrating the theoretical results. 
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1. Introduction. In recent decades methods for approximation of a continuous 
function y on the sphere § 2 := {x = (x\, X 2 , Xs) T el 3 : x\ + x 2 + x 2 = l} by means 
of polynomials have been discussed by many authors (see, for example, jj] [31] 39, 
HO]). Often the underlying motivation has been the need to approximate geophysical 
quantities. For example, such a task appears in the satellite gravity gradiometry 
problem (SGG-problcm) [7], p. 120, 262, [28], in which the task is to find a spherical 
harmonic representation of Earth’s gravitational potential from satellite observations. 
The present study was motivated by this example. We shall return to it several times 
throughout the paper. 

The mathematical problem considered in this paper is to find a polynomial 
approximation to y £ C(S 2 ), given noisy data values y e (xi) at points x, £ S 2 , 
i = 1, ...TV, using a least-squares strategy developed in [l]. (In the SGG applica¬ 
tion the sphere in question is determined by the satellite orbits. The gravitational 
potential at the satellite height is smoother than at earth’s surface, a complicating 
feature for the inverse problem.) We shall assume, in a slight generalization of [T], 
that the point set Xjy :={xi,...,xjv} consists of the points of a cubature rule which 
is exact for all polynomials p £ P 2 m, where Pm is the set of all spherical polynomials 
of degree less than or equal to M, or in other words the restriction to § 2 of the set 
of all polynomials in R 3 of degree less than or equal to M. Thus the point set must 
satisfy 


( 1 . 1 ) 


N 

vp £ p 2 m, y> iP (xQ 

i— 1 


p(x)dw(x), 


where dw(x) denotes area measure on 8 2 , and Wi,i = 1 ,,N are positive cubature 
weights associated with the pointset Xjy. For sufficiently large N one can find in the 
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literature a variety of suitable cubature formulas (see, e.g., mmmm)- Moreover, 
in principle the point sets for such a rule can be generated by selecting from any 
sufficiently dense set of points on the sphere, see EMMS]. 

The strategy is to take the approximant \)m £ Pm to be the minimizer of the 
regularized discrete least-squares problem 



N 


(1.2) y M = argmm l ^Wiipixi)-y £ (xi)) 2 +a'^2'Wi(R M p(.Xi)) 2 , P £ Pm / > 


where y e (x-i) := y(xi) + represent noisy values of a perturbed version y e of the 
original function y calculated at the points of Xn, a is a regularization parameter, 
and Rm : Pm —> Pm is a linear “penalization” operator given by 



(1.3) 



where Pk is the Legendre polynomial of degree k, and in the last step we used DD- 
Here the numbers /3k, k = 1 ,,M are a non-decreasing sequence of positive parame¬ 
ters. With /3o fixed in some appropriate way, the important feature of the parameters 
/3k is their rate of growth. The central task in this paper will be to assign appropriate 
values for the (3k- 

As pointed out in Tj, the expression in m is the most general rotationally 
invariant expression for a linear operator on the space Pm- 

In [T] the point set Xn was taken to be a spherical 2M-design, which simply means 
that HU must hold with equal weights Wi = 4n/N. We gain considerable freedom in 
this paper by allowing general positive weights Wi in DU- The only effective difference 
in the present approximation scheme is that the least-squares problem is slightly 
non-standard because of the appearance of the cubature weights Wi. 

It was observed in numerical experiments in [Tj that a proper choice of the penal¬ 
ization operator Rm together with the regularization parameter a can significantly 
improve the approximation. However, the choice of the model parameters in (11.311 was 
not settled, and still remains an open issue. In our paper we will tackle this crucial 
question by proposing parameter choice strategies (strategies for choosing /3k and a) 
that allow good approximation of noisy smooth functions on the sphere. 

The paper is organized as follows: in the next section we present necessary pre¬ 
liminaries, and give an explicit solution of the regularized least-squares problem.. In 
Section 3 we derive theoretical error bounds for the resulting approximation. Sec¬ 
tions 4 and 5 discuss error bounds and parameter choice strategies. Finally, in the 
last section we present some numerical experiments that test the theoretical results 
from previous sections. 

2. Preliminaries. We introduce a real spherical harmonic basis for Pm, see [251 


{Y kJ : k = 0,1 ,..., M, j = 1,..., 2fe + 1} , 
assumed to be orthonormal with respect to the standard L 2 inner product, 
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Then for p £ P^f an arbitrary spherical polynomial of degree < M there exists a 
unique vector 7 = ( 7 kj) £ ]gd M +i) such that 

M 2fc+l 

(2.1) p(x) = J2 J2 x G §2 ' 

k —0 j =1 

The addition theorem for spherical harmonics (see [22]), which asserts 

2h 1 

(2-2) ^2 Y k>j (x)Y kij ( z) = —-Pfe(x-z), x,z £ S 2 , 

3 =1 ^ 

will play an important role. 

The assumption that a function y on the unit sphere is continuous implies that 
y £ L 2 (S 2 ), and hence that its Fourier coefficients (Y k j.y) L - 2 , s2 , with respect to the 
basis of spherical harmonics are square-summable, i.e. 


00 2 /c+l 


5Z Z | 

k =0 i=l 


L 2 ( S 2 ) 


2 

< OO. 


To measure any additional smoothness of y it is convenient to introduce a Hilbert 
space that is especially tailored to the particular problem, namely 


y£W^ \ 3 : llsll 


2 oo^2H-l |(Y' fcj . ! ^) i2 ( g2 ) 

:= 

k —0 j—1 


W) 


< 00 


where </> is an non-decreasing function such that </>( 0 ) = 0 and P = {/3o, Pi ,..., Pm, ■ ■ ■} 
is the sequence of coefficients appearing in the regularizer (11.31) . In the literature, see, 
e-g., H3, the function goes under the name of smoothness index function. 

In this context the smoothness of y is encoded in </> and p. For example, if the 
smoothness index function (f>(t) and the sequence P = {Pk} increase polynomially 
with t and k such that </>(f) = O [t Vl ), pk = O (k U2 ) ,v\v 2 > 1 / 2 , then the space 
becomes a spherical Sobolev space H 2viV2 (see, e.g., [7], p. 64), and a spherical analog 
of the fundamental lemma due to Sobolev (see [7], Lemma 2.1.5) says that H 2viV2 is 
embedded in the space C^\S 2 ) of functions, which have v continuous derivatives on 
S 2 ,is < 2v\i/ 2 — 1, and are the restrictions to § 2 of functions satisfying the Laplace 
equation in the outer space of § 2 and being regular at infinity. Then Jackson’s theorem 
on the sphere (see (3(Jj, Theorem 3.3) tells us that for y £ W^’^, there holds 


(2.3) 


inf 

P&Pm 


II y .pliers 2 ) — o ) 


v < 2v\ v 2 — 1 . 


On the other hand, if the sequence p = {pk} increases exponentially then for 
polynomially increasing (f> and y £ we have 


inf 

pgp m 


II V alleys 2 ) — 


O ( e -® M ) , 


where q is some positive number that does not depend on M. 
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In the error analysis later in the paper we make use of a linear polynomial ap¬ 
proximation that in a certain precise sense mimics best approximation in the space 
of spherical polynomials of half the degree. The approximation, studied in mmm, 
approximates a function y G C(S 2 ) by VmV € Pm defined by 

m / b \ 2k+1 

(2-4) VmZ/( x ) Y, Ykj (x) (Yfcj, y) Tj 2 (^ 

fc=o ' ' i=i 

= g'* ( u j "tT i Pl(x ■ 

where h is a real-valued function on R + , called a filter function, which satisfies 

(2.5) h(t) g [0,1] Vt e R+, fc(*) = | J’ * | ^ ] ; 

It is shown in [35] that for suitable choices of the filter h (including any filter in 
C 3 (R + ), or the unique C 1 quadratic spline with breakpoints at 1/2, 3/4 and 1 that 
satisfies (ESI), the norm of the operator Vm as an operator from Pm to C(§ 2 ) is 
bounded independently of M. Since, as is easily seen, Vm reproduces polynomials of 
degree less than or equal to M/2, it follows in the usual way that 

II V - v My\\ c ^) < c p j^ lb -Pllc(s 2 ). 

where [•] denotes the floor function. (In this paper c is a generic constant, which 
may take different values at different occurrences.) In view of (12.31) . for polynomially 
increasing <f>, /? and y G we have 

II y - VmvWcw < c [M/2 p < cM~r 

On the other hand, for exponentially increasing /3 and polynomially increasing <fi 
the theory 132] suggests taking h(t) = 1 for t G [0,1] (in which case VmV is just the 
Mth-degree partial sum of the Fourier-Laplace series). Then for y G W ^ there holds 

\\y- VmvWc^ < ||j/ — pII^s 2 ) < cVMe~ qM . 

3. Weighted regularized least-squares problem and its solution. The 

penalization operator El can equivalently be written, using the addition theorem 
(12.21) and (12.11) . as 

M 2fc+l 

(3.i) rmp { x ) = y ] Pk y ' i/c,i( x ) Q^kjiP) 

k —o j =1 
M 2fc+l 

= 7 k,jY kJ (x), 

k —0 j— 1 

allowing us to write the minimization problem as one of linear algebra. For the noisy 
function y e defined on § 2 , let y e := y e (X/v) be the column vector 


y e = [2/ c ( x i), •••,y e ( x Jv)] T GR n 
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and let Ym := Ym(V/v) G R( m+1 ) 2x7V be the matrix of spherical harmonics evaluated 
at the points of X jy. Using this notation we can reduce the minimization problem in 
m to the following discrete minimization problem: 


(3.2) 


min 

T £]R(M + 1)2 


W 1/2 Y m T 7- W 1/2 y e 


2 

2 


+ a 


W 1 / 2 R m T 7 


2 

2 


a > 0 , 


where ||-|| 2 is the standard Euclidean vector norm, Rm := Rm(-^at) = BmYm G 
^(m+i) 2 xAT^ - g a positive diagonal matrix defined by 


(3.3) B m := diag(/3o,^ 1 ,/3 1 ,^ 1) —,Pm) G M (m + 1)2x(m + 1)2 ) 

V -V.-' --v-' 

3 2M+1 

and W is a diagonal matrix of cubature weights 

W := diag(wi,..., w N ) G R NxN . 

The solution of (11.21) can be found from the following system of linear equations 

(3.4) (Y m WY m T + aB M Y M WYM T BM)7 = Y M Wy e . 

Theorem 3.1. Assume y e G C^S 2 ). Let M > 0 be given, and let II.Ill hold true 
for the set of points Xn■ Then j3-4\ ) has the unique solution 7 = (7 k,j) G r( m+1 ) , 

1 N 

(3-5) 7 kj = x + n 2 X] m Y k,j(xi)y e {xi), 

i= 1 

and the minimizer of WM is given by 

M 2/e+l v / \ N 

(3.6) j/m(x) = T^ M y e ft) : = ^ ^ ^ * 2 ^ w i Y ktj {x i )y e (x l ) 

k=0 j=l UPk i= 1 

M or. 1 1 I A 

= £ — ttt? E • x *‘( x T 

/c—0 i= 1 

Proof On using m we have 

N 

^ or, 'y k.j (x v ) U,-(x/ j — ftk,j ; Y K;?t )^ 2 ^g 2 ^ — dk,K,dj,n 

i=l 

where fc, k = 0,..., M, j = 1,..., 2fc + 1, i = 1, ...,2 k + 1. Thus YmWYm T is the 
identity matrix. Since Bm and W are diagonal matrices, the solution of (13.41) is given 
by (13.51) and from (12.11) we obtain (13.61) . □ 

Remark 3.1. Note that one can also employ fast iterative algorithms for scat¬ 
tered least squares [ 12 ] to find the minimizer & Moreover, the evaluation of the 
coefficients (13.51) can be realized with fast spherical Fourier transform presented in 
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4. Error bounds. In this section we estimate the uniform error of approxima¬ 
tion of y by yM, see (13.61) . It is convenient here to regard y e as a continuous function 
on § 2 , constructed by some interpolation process from its values on the discrete set 
Xm- The operator M defined in (13.61) can then be considered as an operator on 
the space C(§ 2 ). Since yM = M y e it is clear that 

y-VM = y- T^ M V M y + T^ M (V M y -y + y- y e ), 


and hence 

(4.1) \\y — 2/m|Ic(S 2 ) < V~ T a,M V My 


t: 


a,M 


C(S 2 ) 


C(S 2 ) 

(\\y ~ ^ / Af2/|lc , (s 2 ) T || y — 2/ e |lc(s 2 )) ’ 


where 


rpft 

± a,M 


is the norm of the operator Tjf M : C( S 2 ) -A C(S 2 ). 
C( S 2 ) ’ 


Let e = [ei, £ 2 , •••, £ a] £ M. N , and ||e||^ = max|ej|. It is natural to assume, and 
from now on we shall do so, that \\y — 2/ e |lo(s 2 ) = II e lloo■ This means that we adopt 
the deterministic noise model, which allows the worst noise level at any point of Xfj. 
Then it is also natural to assume that M is large enough to ensure ||y — VmV^c(% 2 ) — 
Hell^, since otherwise data noise is dominated by the approximation error and no 
regularization is required. Then the bound (ED can be reduced to 


(4-2) ||y yM|| c( s 2 ) < 


V - T a,M V MV 


C( S 2 ) 




,M 


We will call the first term of the right-hand side in (14.21) the regularization error 
and the second the noise propagation error. 

The noise propagation error can be quantified by the following result for the norm 
of M i which is a consequence of (13.61) . 

Theorem 4.1. Under the conditions of Theorem 1 3. 1\ 


(4.3) 


L a,M 


N 

= max ) Wi 

C(S 2 ) XGS 2 “ 
i—l 


M 

E 

k=0 


2k + 1 


4tt( 1 + a/3 2 ) 


Pk(x ■ x. 


N M 2k + 1 

< maxy Wi V . -|Pfc(x • Xj)|. 

^ 4 tt ( 1 + aPl)' 


Theorem 14.11 reduces to Proposition 5.1 in pQ on setting u\ = 47r/iV, but note 
that the result as stated in [T| corresponds to the upper bound in (14.31) . and so is not 
correctly stated. 


Now we are going to bound the regularization error 
start with the following decomposition 


y - T a,M V My 


. We 

C( S 2 ) 


(4.4) 


y - T a,M VM y = y~ T ° ,m v mv + (t 0 ,m - Tf iM )F M y, 
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where T 0j m is the so-called hyperinterpolation operator fM] . 

M 2k+l N 

(4-5) T 0M g (x) = X] Yk >3 ( x ) Y mYk,j(xi)g(xi) 

k—0 j =1 i—1 

M , N 

= Y — Y WiPk ^' x *)£K x *)- 

k—0 i=1 


From (14.51) and CD it immediately follows that for any p £ Pm we have Tq^mP = 
p. Therefore, Tq^mVmV = YmV■ In view of this property and the decomposition (14.41) 
we can derive a bound for the regularization error 


(4.6) y - T^mVmV < II V ~ YmvWc^ + ( T 0 ,m ~ T^ M )V M y 


C( S 2 ) 


< Iklloo + {Tq,m — T a m )VmU 


C(S 2 


An estimate of the term 
theorem. 


{To,m ~ T^ m )Vm 


C{ S 2 ) 


in (14.61) is given by the following 


Theorem 4.2. Assume that the smoothness index function 4> is such that the 
function t —> t/<f>(t) is monotone. Then for y £ there holds 


(4.7) 


{T 0 ,m - T^ M )V M y 


< cMf)(a) \\y\\ w <,,p , 

C (S' 4 ) 


where 4>(a) = 4>(a) if t/<fi(t) is non-decreasing, and <j>{a) = a if t/(f>(t) is non¬ 
increasing. 

Proof. In view of (13.61) . (14.51) and (12.41) . together with the fact that the cubature 
formula in CD is exact for p £ P 2 m, we may write 


(To,m ~ T a M ) V My 


C(S 2 ) 


M 2/c+l 


C( S 2 ) 


Y. Y. * {Yk,j,V M y) L 2^) 

k —0 3 = 1 Pk 

M 2fc+l / j \ n2 

YY hYk d^ x + Y ( y k,j,y) L 2 (S 2) 

k =0 7=1 v k 


where in the last step we used (Yfcj, Vm2/)l 2 (s 2 ) = h(k/M)(Ykj,y) l 2 {% 2 )- Now using 
the Nikolskii inequality (see, e.g., [22], Proposition 2.5) and also h(k/M) < 1, we 
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obtain 


\\(Tq,M Ta, M )VMy\\c(§ 2 )< cM 


M 2/c+l 


c=0 j= 1 
M 2fc+l 


^ ^ h [m) Yk - J 1 +% 2 k 

k= 0 7=1 V 7 


L 2 (§ 2 ) 


E 5 > £ i 


< cM 


^k =0 j=l 
M 2fc+l ✓ /q2 

V V i 

“ “ V1 + a/3| 

A;—0 j=l v K 




1/2 


a# 


v)l 2 (s 2 ) 


( 


<m- 2 ) 


0 / 'k,j,y) L 2 (§ 2 ) 


2 \ 1/2 


W") 


< cM sup 
m 6[0^ o "' 


a 


a + m 


0(w) 


II//II W*’’/ 3 — cM0(a) ||y|| 


where the last inequality follows from 113 . Proposition 2.7. □ 

It is instructive to note that if, for example, 0(t) = t v , then the function 0 defined 
in Theorem 14.21 is given by 


4 >( a ) = 


a, v > 1, 
a", 0 < v < 1. 


Thus the error bound in the theorem does not improve if 0(t) grows faster than t. 

5. Parameter choice strategies. In this section we will be concerned with the 
choice of the design parameters for the least-squares approximation dm, namely the 
regularization parameter a and the penalization parameters f3k- In the first subsec¬ 
tion we discuss an a priori choice for the penalization parameters /3k- In the next 
subsection we consider an adaptive strategy for choosing the regularization parame¬ 
ter a. In the third subsection we present an a posteriori choice for the penalization 
parameters /3k- 

The choice of parameters is motivated by the error bound (14.211 for y — yM- From 
mm and (14.71) it follows that the bound (14.21) can be reduced to the following: 

(5.1) ||y -Vm\\< IklL + cM0(a) \\y\\ w4 ,,f> H 

< cM0(a) \\y\\ W 4,,f> TcHell^ 

5.1. A priori choice of the penalization parameters. For definiteness, we 
assume in this subsection that 0(i) = t, which means that 0 has the highest order 
in a, namely 0(a) = a. The error bound (15.11) now provides useful guidance in the 
choice of the regularization parameters /3k- If /3o is considered to be fixed, and we 
increase the rate of growth of the /3k, then the first term on the right-hand side of 
the last line of JMD will increase, while from (14.31) the second term has an upper 
bound that decreases with increasing rate of growth of the /3k- Even more can be 
said: for the first term to be finite the norm of y must be finite, which imposes 

the constraint 


■2 Hell 


1 a,M 


rpf: 




C{ S 2 ) 




oo 2fc+l 


yi yi Pk { Y k,j,y) L 2(&) 

k =0 j—1 


(5.2) 


< oo. 
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To see what this condition means in a particular application, we consider the SGG- 
problem mentioned in the Introduction. In this problem y is the second order radial 
derivative of the gravitational potential measured pointwise at the orbital sphere of a 
satellite. It can be shown Emm that after a proper normalization of this sphere 
to 8 2 we have 


(5.3) (^fe,j>2/)z,2(s2) ~ a k9k,ji 

where a*, = ( 77 ) (fc+!)(fc+ 2 ) , ^ rac ji us 0 f the orbital sphere, R is the radius 

of the surface of the Earth considered as a sphere, and is some (unknown) 

sequence of scaled Fourier coefficients of the gravitational potential g measured at 
the surface of the Earth. It is well-known (see, e.g., 1330 ) that in the scale of the 
spherical Sobolev spaces {H s } mentioned above the Earth’s gravitational potential has 
a smoothness index s = 3/2 at least, which means that the sequence {gk,j} should 
satisfy the requirement 


00 2 fc+l 

+ 1 / 2 ) 3 9k,j < oo- 

fc =0 j =1 

In view of the last requirement, the condition (15.21) is satisfied by the choice 

(5.4) /3k = a / 7 1/2 (fc + l/2) 3/4 , A; = 0,1,.... 

Of course the condition (15.21) will also be satisfied if the /3k increase more slowly, but 
at the likely expense of a larger second term in the error bound (15.11) . 

Since R < p, it is clear that the (3k given by (15.41) increase exponentially. This 
is natural in view of the exponential decrease of the Fourier coefficients (15.31) of the 
approximated function, which implies that the exact function as measured at the 
satellite height is very smooth, even analytic. The regularization scheme with 
weights (15.41) will penalize the presence of oscillating coefficients with large indexes in 
the approximant M y e - In the last section we illustrate a good performance of the 
scheme (EH) with these penalization weights. 

5.2. Regularization parameter choice strategy. For regularization of our 
problem we will implement an adaptive regularization parameter choice strategy 
known as the balancing principle (see, e.g., [T7] HH ;29] and references therein). 
In this method the regularization parameter a is selected from some finite set, say 
A l := {oj = q l a 0 , *=1,2,..., L}, with q £ (0,1) and L large enough. 

Applying the balancing principle to our problem we start with the smallest pa¬ 
rameter cxl and increase stepwise a^_ 1 = cti/q,i = L,L — 1,..., until a* := a z is the 
parameter for which 


Ta z ,M'y t T o 2 + 1 ,m!/ 


C{ S 2 ) 


> ae e 


T‘ 


a z +i,M 


C(§ 2 ) 


for the first time. Here ae is a design parameter. In all our numerical tests with 
the balancing principle (BP) reported below in Section 6 , the value of ae is fixed as 
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se = 0.002 and is data independent, while the value of the regularization parameter 
a* chosen according to BP varies with data. 

Note that for choosing a* we need only the knowledge of (13.61) and an upper 


bound of 


T: 


ot,M 


C{ § 2 ) 


given by (14.31) . 


In the Section 6 we will present a numerical test showing a good reconstruction 
of the function on the sphere from noisy observations with the above a posteriori 
regularization parameter. It is instructive to see that in all tests BP performs at the 
level of the ideal parameter choice a G A^. 


5.3. A posteriori choice of the penalization weights. We start with the 
observation that the space Pm of spherical polynomials p is a reproducing kernel 
Hilbert space (RKHS) W. By the Riesz representation theorem, to every RHKS H 
there corresponds a unique symmetric positive definite function K : § 2 x S 2 -> I, 
called the reproducing kernel of H = Hk, that has the following reproducing property: 
p(x) = (p(-), K(-,x)) nK . A comprehensive theory of RKHSs can be found in [2|. 

It is easy to check that the kernel 


M 2k+l 

(5.5) K{x, z) = J2&k 2 £ y ‘j( x )^J'( z )> x,z <E § 2 

fc=o i=i 

has the above mentioned reproducing property if the inner product in Pm is defined 
as follows 


M 2fc+l 

(Xkjlf)^ 2 (§ 2 ) {Yk,j,g) L 2 (g 2 ) • 

k —0 j —1 


Indeed, for p € Pm we find 


M 2fc+l 

(p(‘)i ( x ’ '))hk L 2 ^) 0^k,j, K(x, ■))l2(§2) 

k—0 j —1 

M 2fc+l 

= £ (^j,p) i ,2 (S 2 ) ^ 2 n-j(x)=Kx). 

k =0 j —1 

In this RKHS setting the spherical polynomial j/m = Dm (N - , K , a) defined by 

0 also can be seen, using the addition theorem and 0, as the minimizer of the 
following quadratic functional 

N 

(5.6) T a (N,K;p) = ^ Wi(p(xi) - y e (*i)) 2 + a \\p\\ 2 Hk , p £ Pm, 

i=l 

which makes (15.51) a natural way of defining the reproducing kernel in this context. 

At this point the problem of the choice of the penalization weights {fik} is trans¬ 
formed into that of selecting a kernel K from the set /C of kernels of the form (15.51) . 

In the literature there are several methods for choosing a kernel from the available 
set of kernels (see, e.g., [221 27] , and references therein). For example, in [22] the 
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authors suggest selecting a kernel by minimizing the value of the functional (15.61) 
evaluated at its minimizer Dm■ In the present context, according to [52], the kernel 
K = K * of choice is given as 


K * = argmin {T a (N,K;y M (N,K,a)), K € 1C} . 

Note that such K * depends on the value of the regularization parameter a. Therefore, 
the approach 122] can be realized only for an a priori known a. However, in practice 
we are not provided with this knowledge, and have to use a posteriori regularization 
parameter choice strategies (for example, the balancing principle described in Sub¬ 
section 4.1). Thus, in practice we are dealing with kernel dependent regularization 
parameter a = a(K). 

This situation has been discussed in m- In the present context the kernel choice 
suggested in [27 can be written as follows 


(5.7) K+ = argmin {T a ^ K ^ {N , K;y M (N , K, a(K))), K € Kj . 

The existence of such K + has been proved in m under rather general assumptions on 
the set of admissible kernels 1C and regularization parameter choice strategy a = a(K). 

From a practical point of view, it is a challenging issue to use the strategy m 
in our case because one has to minimize a function depending onM + 1 unknown 
penalization weights /3fc. Therefore, it is natural to reduce the complexity of the model 
before applying the strategy from m- 

For example, one may parametrize {(3k} as follows: /3 2 = e Al ^ fc+1 ^(fc+l) A2 , Ai, A 2 > 
0. In other words, in (15.71) the set of kernels K, consists of the functions 

M 2k+l 

(5.8) K(t,r) = J2e~ Xl{k+1) (k + l)- x * ^ Y kjj (t)Y kd (r), t,r € S 2 . 

fc=0 j=1 

Then the kernel K + can be found by minimizing a function of two variables A 1; A 2 . 
In the last section we will illustrate such a reduced approach by a numerical test 
showing good performance of the scheme m with a posteriori chosen penalization 
weights. 

6. Numerical examples. In this section we present some numerical experi¬ 
ments to verify the analysis from the previous sections. Note that we work not with 
real data but with artificially generated ones. In all our experiments we follow mm 
and assume that the set of points Xn is the set of Gauss-Legendre points, for which 
the positive quadrature weights are known analytically. The number of points in this 
case is N = 2(M + l) 2 , and the corresponding cubature formula (11.11) is indeed exact 
for all spherical polynomials of degree 2 M. In all our experiments M = 30. 

Note that in real applications the spherical polynomials of much higher degree 
are used [25j. Moreover, the Gauss-Legendre points are known to have the drawback 
of having too many points concentrated at the poles, making it not suitable for real 
satellite data. In our experiments below we use the Gauss-Legendre points and poly¬ 
nomials of modest degree only for illustration purposes and as a proof of concept. At 
the same time, we note that even for the case M = 30 the corresponding discrete 
problem is rather ill-conditioned and, thus, should be treated with a regularization 
(see Figure I7TT1 and the discussion below). 
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We start with an experiment illustrating that a proper choice of the penalization 
weights /3q , Pm is crucial for the approximation of functions on the sphere. Consider 
again the SGG-problem corresponding to (15.31) . Note that for k = 1, 2,30 the values 

ak = ("f) i fc+1 ^ fc+2 i in (15.31) are increasing, and so, they do not exhibit a typical 
behavior of the singular values of the compact operators. This effect is well-known 
(see, e.g, [7], Fig. 4.2.3, p. 280). 

Therefore, to mimic the SGG-problem for M = 30 one usually omits the factor 
(fc+ i )(fc + 2 ) ( gee ^ e _g ^ [ 4 j) i n thi S case the decay character of the coefficients ak in 
(15.31) can be modeled, for example, as 

a fe = (1.2)- fe , fc = 0,1,..., M. 

We conduct our first experiment in the following way. First we generate a spherical 
function 


M 2k+l 1 

y = 2/W = ^(i.2)“ fe Y 9 k,j—Y k ,j{x), x e s 2 , 

fc=o j=i P 

where gkj = (fc + 1/2 )~ 3 ^ 2 Xk.j, k = 0 ,..., M, j = 1 , ...2k + 1 , and x k j are random 
numbers uniformly distributed on [0,1]. The blurred spherical function y e is simulated 
by adding a random point-wise noise to the values of the initial function y at the point 
set Xn- The simulated noise values are given as the components of a random vector 
0.05e/ Hell^, where e = [ei, e 2 ,..., cat], and et are uniformly distributed on [—1,1], To 
mimic the SGG-problem we reconstruct the vector g = (gk,j) by g a ’ M = (r//') Vf ), 
where g^’J 1 = a^jk.j, and 7 k,j are given by (13.51) . 

To assess the obtained results and compare the performance of the considered 
schemes we measure the relative error 


lldll 2 ‘ 

The results are displayed in Figure IG.ll where along the vertical axis we plot the 
relative errors in solving the problem with one of 50 simulated data. The relative 
errors are plotted in ascending order for each of four methods: a straightforward 
least-squares fit to noisy data without any regularization, the regularization with the 
penalization weights (15.41) and a chosen according to the balancing principle (BP) 
from Ago = {cq = 8 • q l i = 1 , 2 ,..., 60} , q = 0 . 8 , the regularization with default 
penalization weights /3^ = 1,A; = 0,1,..., M, and the best a £ Ago, the regularization 
with the penalization weights (15.41) and the best a £ A 60 . Thus, in the latter two cases 
the choice of the regularization parameter a for both schemes was made to achieve 
the best possible performance of each method. As it can be seen from Figure IfTTl the 
balancing principle (BP) performs at the level of the ideal parameter choice strategy. 

From Figure 16.11 one can also conclude that the proper choice of the penaliza¬ 
tion weights according to the proposed a priori recipe can significantly improve the 
accuracy of the reconstruction. Moreover, Figure 16.11 shows that a straightforward 
least-squares fit to noisy data without regularization leads to the relative error that 
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Fig. 6.1: Numerical illustration. The figure presents relative errors for 50 simulations 
of the data. The errors are plotted in ascending order for each of the discussed 
methods. Note that two bottom curves corresponding to penalization according to 
(15.41) nearly overlap. 


is about 2-3 times larger than that after a regularization. This confirms that in the 
considered experiment we are really dealing with a rather ill-conditioned problem. 

In our second experiment we again confirm that the balancing principle gives a 
value of the regularization parameter a* that is competitive with the best value man¬ 
ually found in T] . We choose the regularization parameter from the same geometric 
sequence Ago and use the same value of the design parameter ae = 0.002 in BP. 

Similarly to [Tj, as a test function y we take the sum of the Franke function y\ 
modified by Renka |23] (p.146) and a function y cap [Hj . namely y = yi + y c ap with 


y 1 (x 1 ,x 2 ,xs) = 

+ 0 75 e -( 9a; i+ 1 ) 2 /49-(9x2 + l)/49-(9x3 + l)/10 

+ 0.5e _(9a:i_7)2/4_(9::C2_3)2/4_(9a:3_5)2/4 

- 0 . 2e -( 9 - 1 - 4 ) 2 -(^-7) 2 -(9x3-5) 2 ) ( XuX2}X3 ) e §2, 


and 


( 6 . 2 ) 


2/cap (X) 


2 cos (n arccos(x c • x)), x c • x > cos(0.5), 
0, otherwise, 


where x c = ^,^=J and (•) defines the dot product of two vectors. The 

function y was then contaminated by noise, taking for the noise e(x) at each x £ Xjy an 
independent sample of a normal random variable with mean 0 and standard deviation 
cr = 0.5. 
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(b) y with N( 0, 0.25) noise 


(a) V 



(c) T a *,My e 


Fig. 6.2: Franke function recovery 


Figure [fOa] illustrates the function y, while Figure [6~.2bl shows the blurred function 
y e (x) =y(x) + e(x). 

For the reconstruction, following [Tj we choose a Laplace-Beltrami penalization 
operator that corresponds to the matrix 


B m := diag(0, 4 , 4 , 4 , (M(M + l )) 2 , {M(M+1)) 2 ) £ r(m+i) 2 x(m+i)^ 

3 2M+1 

Figure IQcl illustrates the reconstructed function T f M y e . The regularization 
parameter a* was obtained according to the balancing principle described above. We 
found automatically the regularization parameter a* = 1.42 • 10~ 4 which agrees well 
with the value 10 -4 from [I] obtained manually. 

In our last experiment we will illustrate an application of the a posteriori rule 
(EH for choosing the penalization weights. As a test function y e we again consider 
the blurred function from the previous example, where we used the a priori chosen 
penalization weights /3k = k(k + 1) corresponding to Laplace-Beltrami operator. Now 
we are going to estimate the penalization weights using the a posteriori strategy 
described in Subsection 4.3. 
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Fig. 6.3: Numerical illustration. The figure presents relative errors for 50 simulations 
of the data. The errors are plotted in ascending order for each of the discussed 
methods. 


Recall that we are looking for the minimizer (15.711 among the set of admissible 
kernels K, consisting of the functions (15.811 . This approach allows us to take into 
account an exponential, as well as a polynomial growth of /3k- 

To find an approximate minimizer of (15.711 we have implemented the Random 
Search method [19] over the set of parameters (Ai,A 2 ) € [0,5] x [0,5]. The method 
was implemented 10 times, and in each implementation 10 random steps have been 
performed. Then the mean values of the parameters Ai,A 2 appearing after each 
implementation of the Random Search method have been taken as an approximate 
minimum point. As the result, the values Ai = 0.32, A 2 = 1.9 have been obtained. 

Figure [6~3l displays the relative errors in solving the problem (16.11) . (16.21) with one 
of 50 simulated noisy data, for each of two methods: regularization with the penal¬ 
ization weights Pk = k(k + 1), and regularization with a posteriori chosen weights. 

From Figure 15731 we see that the choice of the penalization weights according to 
the proposed a posteriori choice rule can improve the accuracy of the reconstruction. 
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